Solution of the isothermal sphere in hydrostatic Equilibrium

This is derivation of the gravitationally bound sphere (1D) in hydrostatic equilibrium for the isothermal case. The goal is to get a density distribution of the sphere in 1D depending on the radius $r$

First, we begin with the hydrostatic equilibrium:

\begin{equation} dP = -\rho g dr \end{equation}

of course, $\rho $ and $g$ are dependent of radius and the pressure $P$ is a function of density $\rho$

$$dP(\rho) = -\rho(r) g(r) dr$$

Now we use the equation of state for an isothermal gas:

$$ P = \rho R_s T \\ dP = d\rho R_s T $$

We insert this into the above equation of the hydrostatic equilibrium.

Other useful relations are:

\begin{equation} d\rho R_s T = \rho(r) g(r) dr \\ g=G\frac{M}{r^2} \\ dM=4\pi \rho(r)r^2dr, \qquad \frac{dM}{dr}=4\pi \rho(r) r^2 \end{equation}

We can write $d\rho / \rho$ as a logarithmic derivative:

$$\frac{d\rho}{\rho} = -\frac{GM(r)}{R_sT}\frac{dr}{r^2} \\ d\ln \rho = - \frac{GM(r)}{R_sT}\frac{dr}{r^2}$$

And some more sorting of the terms in the equation:

$$r^2\frac{d\ln \rho}{dr} = -\frac{G}{R_sT}M(r)$$

Doing a derivative, $d/dr$ and inserting what we have for $dM/dr$ we get:

$$ \frac{d}{dr}\left( r^2 \frac{d\ln\rho}{dr} \right)= -\frac{G}{R_s T} \frac{dM(r)}{dr} \\ \frac{d}{dr}\left( r^2 \frac{d\ln\rho}{dr} \right)= -\frac{4\pi G}{R_s T} \rho r^2 $$

This is the differential equation for the isothermal sphere. Now to the solution wwhich has to be obtained numerically.

To solve this differential equation, we have to do some algebra

We need to seperate out the first order derivatives and the second order derivatives. To do this we apply the product rule, which means that the differential operator acts on all the terms seperately. If we factor this out, we get a non-linear differential equation of second order. Differentiation $d/dr$ is denoted with the ' sign for shorthand notation.

$$ \frac{d}{dr}\left( \frac{r^2}{\rho} \frac{d\rho}{dr}\right) = -\frac{4\pi G}{R_s T} \rho r^2 \\ \frac{d}{dr}\left( \frac{r^2}{\rho}\right) \frac{d\rho}{dr}+ \frac{r^2}{\rho}\frac{d^2\rho}{dr^2} = -\frac{4\pi G}{R_s T}\rho r^2 \\ \frac{2r}{\rho}\frac{d\rho}{dr}-\frac{r^2}{\rho^2}\left(\frac{d\rho}{dr}\right)^2 + \frac{r^2}{\rho}\frac{d^2\rho}{dr^2}= -\frac{4\pi G}{R_s T}\rho r^2 \\ \frac{2r}{\rho}\rho' - \frac{r^2}{\rho^2}\rho'^2 + \frac{r^2}{\rho}\rho''= -\frac{4\pi G}{R_s T} \rho r^2 \\ \frac{2r}{\rho}\rho' - \frac{r^2}{\rho^2}\rho'^2 + \frac{r^2}{\rho}\rho''+\frac{4\pi G}{R_s T} \rho r^2 = 0 \\ \rho'' - \frac{\rho'^2}{\rho} + 2\frac{\rho'}{r} + \frac{4\pi G}{R_s T}\rho^2 = 0 $$

Substitutions to make the equations dimensionless in order to be integrated numerically:

$$ y=\frac{\rho}{\rho_c}, \qquad x=r\sqrt{\frac{4\pi G\rho_c}{R_s T}}, \qquad dr=\frac{dx}{\sqrt{\frac{4\pi G \rho_c}{R_s T}}} $$

The derivative $d/dr$ changes to $\frac{d}{dx}\sqrt{\frac{4\pi G\rho_c}{R_s T}}$ and eliminates the prefactor $\frac{4\pi G}{R_s T}$ in the equation. We then have the dimensionless differential equation:

$$y'' - \frac{y'^2}{y} + 2\frac{y'}{x} + y^2 = 0$$

To get this into a system of two first order differential equations, we do the following:

$$ \frac{dy}{dx}= y' \\ \frac{dy'}{dx}=\frac{y'^2}{y} - 2\frac{y'}{x} - y^2 $$

We can use this to integrate the density. SciPy provides a function odeint which takes a vector of derivatives as input and returns the value.



In [22]:

    
%matplotlib inline
plt.rcParams['figure.figsize']=(9.0,6.0)
plt.rcParams['font.size']=16
plt.rcParams['font.family']='serif'



In [23]:

    
import numpy as np
import scipy.constants as C
from scipy.integrate import odeint
import matplotlib.pyplot as plt

N=10000 # number of points on R
T0=200. #arbitraty temperature
rho_c = 1e-10
c=np.sqrt((4*C.pi*C.G*rho_c)/(C.R*T0))

#derivative function which is needed for the integration
def deriv(y,r):
    return np.array([y[1], y[1]**2/y[0] - 2*y[1]/r - y[0]**2 ])

yinit = np.array([1., 0.])       # initial value vector
radius = np.linspace(0.1,100.,N) # radius r=0. cannot be used, as it is divided by r
y = odeint(deriv, yinit, radius) # calling the integration function

fig1= plt.figure(1)
sub1 = fig1.add_subplot(111)
sub1.plot(radius, y[:,0])

sub1.set_xlabel("radius")
sub1.set_ylabel("density $\\rho / \\rho_c$")
sub1.set_yscale('log')
sub1.set_xscale('log')
sub1.set_ylim(1e-2,1.5)
sub1.set_xlim(1e-1,1e2)
sub1.set_title("Logarithmic plot of the Sphere")









    Out[23]:





<matplotlib.text.Text at 0x7fda1e78db00>



In [24]:

    
fig2 = plt.figure(2)

sub2 = fig2.add_subplot(111)
sub2.plot(radius, y[:,0])
sub2.set_xlim(0,20)
sub2.set_title("Linear plot")
sub2.set_xlabel("radius")
sub2.set_ylabel("density $\\rho/\\rho_c$")









    Out[24]:





<matplotlib.text.Text at 0x7fda1eb89160>



In [24]: